# ---
# Fudenberg, Gao, & Liang:
# "How Flexible is that Functional Form? Quantifying the Restrictiveness of Theories"
# 
# Application 3: Microfinance Takeup

# "NetNonlinearComp.R"

# Research assistant: Stephanie Nam
# Date: 09/24/2023

# NOTE: This R script computes the completeness of the partially linear model

#- --

library(paramtest)
library(dplyr)
library(expm)

load(file = "NetNonlinearData.RData")
sink('NetNonlinearComp.log')

# Initialize parameter ranges
theta0_range <- seq(-25, -10, by = 0.1)
#theta0_range <- seq(-3, 0, by = 1)
theta0_range
theta1_range <- seq(0.01, 0.3, by = 0.01)
#theta1_range <- seq(0.01, 0.02, by = 0.01)
theta1_range


# Compute powers of adjacency matrices

A1 <- networks
A2 <- lapply(networks, function(A) A%^%2)
A3 <- lapply(networks, function(A) A%^%3)
A4 <- lapply(networks, function(A) A%^%4)
A5 <- lapply(networks, function(A) A%^%5)
A6 <- lapply(networks, function(A) A%^%6)

# Define functions

calc_H <- function (theta1, A1, A2, A3, A4, A5) {
  H <- theta1 * A1 + theta1^2 * A2 + theta1^3 * A3 + theta1^4 * A4 + theta1^5 * A5 # + theta1^6 * A6 
  return(H)
}

calc_x_j <- function (theta0, H, non, lead) {
  x_j <- c()
  for (i in 1:length(non)) {
    x_j <- append(x_j, theta0 + sum(H[lead, non[i]]))
  }
  return (x_j)
}

calc_p_ij <- function(x_ij) {
  p_j <- 1 / (1 + exp(-x_ij))
  return (p_j)
}

calc_f_i <- function(t_0, t_1) {
  H_i <- mapply(calc_H, theta1 = t_1, A1, A2, A3, A4, A5)
  x_ij <- mapply(calc_x_j, theta0 = t_0, H = H_i, non = nonleaders, lead = leaders)
  f_i <- c()
  for (i in 1:43) {
    f_i <- append(f_i, 
                  sum(mapply(calc_p_ij, x_ij = x_ij)[[i]]) / length(nonleaders[[i]])
    )
  }
  return (f_i)
}


#########################################################################
############################# Completeness ##############################
#########################################################################


print("results")

# MSE between Naive Model and True Take Up Rate
nm_error <- 1/43 * sum((mean(mf_rate) - mf_rate)^2)
nm_error


MSE_lm_func <- function(iter, t_0, t_1) {
  fm <- lm(mf_rate ~ calc_f_i(t_0, t_1))
  e_theta <- deviance(fm) / 43
  return(e_theta)
}


MSE_iter <- grid_search(MSE_lm_func, params = list(t_0 = theta0_range, t_1 = theta1_range))

MSE_min <- data.frame()
MSE_min <- rbind(MSE_min,
                 cbind(MSE_iter[[2]][which.min(unlist(MSE_iter[[1]])),],
                       error = MSE_iter[[1]][[which.min(unlist(MSE_iter[[1]]))]]))

print("range of minimzers")
range_t0 <- c(min(MSE_min$t_0),max(MSE_min$t_0))
range_t0 
range_t1 <- c(min(MSE_min$t_1),max(MSE_min$t_1))
range_t1

print("error of nonlinear model")
err_nl <- MSE_min$error
err_nl

reg_opt <- lm(mf_rate ~ calc_f_i(MSE_min$t_0, MSE_min$t_1))
resid_nl <- reg_opt$residuals
var_err_nl <- 42/43 * var(resid_nl^2)

print("error of naive model")
err_naive <- var(mf_rate)*42/43
err_naive
resid_naive <- mf_rate - mean(mf_rate)
var_err_naive <- var(resid_naive^2)*42/43

print("completeness")
comp_nl <- 1 - err_nl/err_naive # Completeness

cov_err_nl <- 42/43 * cov(resid_nl^2, resid_naive^2)

avar_comp_nl <- 1/ (err_naive^2) * (var_err_nl - 2 * (1-comp_nl) * cov_err_nl + (1-comp_nl)^2 * var_err_naive)
se_comp_nl <- sqrt(avar_comp_nl / 43) # Standard error for completeness

COMP_nonlinear <- data.frame(comp_nl, se_comp_nl)

COMP_nonlinear

save.image(file = "NetNonlinearCompResults.RData")

sink()
